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Betelgeuse, one of the most studied red supergiant stars’’, 
dimmed in the optical range by ~1.2 mag between late 2019 
and early 2020, reaching a historical minimum?” called ‘the 
Great Dimming’. Thanks to enormous observational effort to 
date, two hypotheses remain that can explain the Dimming": 
a decrease in the effective temperature®’ and an enhance- 
ment of the extinction caused by newly produced circumstel- 
lar dust®°’. However, the lack of multiwavelength monitoring 
observations, especially in the mid-infrared, where emission 
from circumstellar dust can be detected, has prevented us 
from closely examining these hypotheses. Here we present 
4.5 yr, 16-band photometry of Betelgeuse between 2017 and 
2021 in the 0.45-13.5 1m wavelength range making use of 
images taken by the Himawari-8" geostationary meteorologi- 
cal satellite. By examining the optical and near-infrared light 
curves, we show that both a decreased effective temperature 
and increased dust extinction may have contributed by almost 
equal amounts to the Great Dimming. Moreover, using the 
mid-infrared light curves, we find that the enhanced circum- 
stellar extinction actually contributed to the Dimming. Thus, 
the Dimming event of Betelgeuse provides us with an oppor- 
tunity to examine the mechanism responsible for the mass 
loss of red supergiants, which affects the fate of massive stars 
as supernovae". 

Himawari-8 is a Japanese geostationary meteorological satellite 
orbiting 35,786km above the equator at 140.7°E (ref. '°). Since 7 
July 2015, Himawari-8 has taken images of the entire disk of the 
Earth once every 10 min using its optical and infrared imager, the 
Advanced Himawari Imager (AHI; Extended Data Fig. 1). The 
Himawari-8 satellite also observes the region of outer space around 
the edge of the Earth's disk during every scan; this motivated us to 
develop a new concept: using meteorological satellites as ‘space tele- 
scopes’ for astronomy (Fig. 1). 

One of the unique aspects of this concept of using Himawari-8 
for astronomy is that it enables us to obtain high-cadence time 
series of mid-infrared images, which are hard to acquire with the 
usual astronomical instruments. Small ground-based telescopes 
can acquire time-series measurements”, but the observations only 
cover wavelengths within the atmospheric window, and they are 
interrupted by the Sun for several months in most cases'*"*. In con- 
trast, astronomical satellites and aircraft-carrying telescopes can 
observe stars in most wavelength ranges'*"’, but the observations 
require higher costs than ground-based telescopes. Survey tele- 
scopes orbiting the Earth for non-astronomical purposes—such 
as the Himawari-8 meteorological satellite—have the potential to 
overcome these problems. In this work, we demonstrate that the 
mystery of the Great Dimming of Betelgeuse can be solved with 
the aid of the time-series observations from Himawari-8, especially 


those in the mid-infrared, with which the amount of dust around 
Betelgeuse can be determined. 

Examining the images taken by the Himawari-8 satellite, we 
found that several bright stars, including Betelgeuse, sometimes 
appear in the images (Extended Data Fig. 2). We therefore mea- 
sured the amount of light coming from Betelgeuse, typically once 
per 1.72d, between January 2017 and June 2021, and we made a 
4.5 yr catalogue of light curves in 16 bands covering 0.45-13.5 jm 
(Extended Data Fig. 3a). Comparing the light curves, that is, the time 
series of spectral energy distributions (SEDs), with a model SED for 
Betelgeuse, we determined the time variation of the main param- 
eters of Betelgeuse, namely the radius R, effective temperature T., 
and extinction A(V) in the optical range. In addition to these three 
parameters, which can be determined from optical or near-infrared 
observation*'*, we also succeeded in determining the dust optical 
depth z,, at 101m using our mid-infrared light curves. With this 
T1o, we Will examine the dust-production history around Betelgeuse. 

Figure 2 shows the main results of this paper, namely, the time 
variations of the stellar parameters R, T., A(V) and 7, together with 
the AB magnitudes for band 3 at 640nm, which is located around 
the central part of the Johnson R band. From the time variations of 
R, Ty and A(V), we found that the Great Dimming of Betelgeuse, 
the drop of ~1.2 mag in the V magnitude’, is likely to be caused by a 
cooling of ~140K in T,,, which results in a dimming of the V mag- 
nitude by ~0.6 mag, and an increase of ~0.6 mag in A(V), which, by 
definition, darken the V magnitude by ~0.6 mag. This result sup- 
ports a previously suggested scenario in which the Dimming was 
caused by a combination of a decrease in T., and an increase in 
A(V) (refs. '*), although several shortcomings in our model SED 
should be kept in mind; some of the shortcomings are examined in 
Methods and Extended Data Fig. 4. 

However, the cause of the enhanced A(V) during the Great 
Dimming is still under debate. Several works have proposed that 
occultation by a newly formed circumstellar dust cloud might be 
a promising cause'’. In contrast to this dust-formation scenario, 
several other works have argued that the dust formation is unneces- 
sary to explain the Dimming. For example, three-dimensional (3D) 
effects’ or molecular opacity’®”, rather than dust extinction, can 
also explain the Dimming. 

To examine whether or not the enhanced A(V) during the 
Dimming is due to the circumstellar dust, we focus on the time vari- 
ation of 7,9. Since 7,) was determined using the mid-infrared light 
curves of bands 12-14 at around 10pm, where O-rich dust emis- 
sions dominate the flux”!”’, it reflects the amount of circumstellar 
dust that emits the infrared light. Therefore, using our infrared 7,, 
measurements, in contrast to optical diagnostics, we can directly 
trace the amount of circumstellar dust. Our measurements in Fig. 2e 
show an enhancement of 7,, during the Great Dimming, though 
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Fig. 1| Betelgeuse observed with the Himawari-8 meteorological satellite near the edge of the Earth's disk. a, False-colour image taken by the 
Himawari-8 satellite on 19 January 2020. Betelgeuse is located at the top right corner of the image. b,e, Zoomed-in images of band 7. d, The response 
functions of Himawari-8; see the basic properties of the Himawari-8 data in Extended Data Fig. 1. e, The observed SED after averaging all observations with 


a model spectrum fitted to it. 


there is a slight time delay between the peaks of A(V) and 7,,. This 
enhancement of z,, indicates that a clump of gas did in fact produce 
dust that shaded the photosphere of Betelgeuse and contributed to 
the Great Dimming. 

Two potential dust reservoirs have been found so far for red 
supergiants. The first is a diffuse and spread-out dust layer located a 
few to several tens of stellar radii from the photosphere””’. In addi- 
tion, instruments with very high spatial resolution have recently 
been used to find circumstellar dust clouds around red supergi- 
ants just a few tenths of a stellar radius above the photosphere**™. 
The distribution and amount of dust in the spread-out layer remain 
mostly unchanged with time”, and thus it is unlikely that dust 
condensation in this layer is related to the large time variation of 
A(V) and 7,, during the Dimming that we found. In contrast, the 
latter, lukewarm region” coexists with the warm chromosphere, the 
temperature of which””’ is sufficiently high for dust sublimation”. 
Therefore, a dust clump that condenses near the photosphere may 
soon be sublimated by chromospheric heating, or, if the radiative 
acceleration is insufficient, it may fall back to the photosphere to 
be sublimated”. Considering these points, we conclude that the 
enhancements of A(V) and 7, during the Dimming may have 
occurred very close to the photosphere. This hypothesis is also sup- 
ported by the fact that the polarized brightness of this circumstellar 
envelope changed substantially around the Dimming”. 

To further examine the dust-condensation process, we focus 
here on the time variation of the amount of molecular gas around 
Betelgeuse, in which the dust condensation occurred, and in par- 
ticular on H,O molecules. In addition to the dust emission seen at 
wavelengths around 10pm, the number of molecules, particularly 
H,0, can also be traced using Himawari-8 photometric data. Since 
H,O molecules in the Earth’s atmosphere absorb starlight” in the 
wavelength ranges where the H,O molecular features from stars 
exist, the observation of extraterrestrial H,O molecules usually 
requires space telescopes'’. One of the missions of the Himawari-8 
satellite is to observe the distribution of water vapour in the Earth’s 
upper atmosphere'’, and thus Himawari-8 has the wavelength 
bands where H,O features exist. These bands, bands 8-10 at around 
6-8 pm, are therefore useful for time-series observations of H,O 
molecules around Betelgeuse. We focus here on the flux variation of 
band 8, where the signal-to-noise ratio of the H,O feature is maxi- 
mized, thanks to the strength of this feature“ and the sensitivity 
of the Himawari-8 satellite (Extended Data Fig. 1). We note that the 


6-8 1m wavelength range is a complicated region in which other 
molecular bands are present, for example, SiO molecules. Thus, we 
should keep in mind the possibility of contamination by these mol- 
ecules, and further spectroscopic observations are required to fully 
disentangle the problem. Nevertheless, band 8 is the least contami- 
nated by other molecules’*”™. 

Figure 3 shows the time variation of the flux excess in band 8, 
defined as the ratio of the observed flux to predicted photospheric 
flux minus 1, which is likely to trace the amount of H,O molecu- 
lar gas around Betelgeuse, maybe in the MOLsphere’***. When gas 
clumps outside the photosphere are present that emit or absorb H,O 
molecular features, the excess is positive or negative, respectively. 
In the first half of the observing period, that is, before September 
2018, to within the error bars the flux excess was almost consistent 
with zero, which indicates that there was no extra detectable H,O 
gas around the photosphere. In contrast, most of the data points 
after September 2018 show negative flux excesses. This indicates the 
emergence of a clump and/or a layer of cool gas containing H,O mol- 
ecules, in which the dust that dimmed Betelgeuse may have formed. 
Considering the clumpy nature of the dust cloud that dimmed 
Betelgeuse’ and the existence of the warm gas clump found from 
ultraviolet observations of the Mg mh and k lines around October 
2019°, the H,O gas that we found is also likely to be a clump. 

Finally, and more interestingly, we tentatively identified a rapid 
transition in the H,O feature from absorption to emission in early 
April 2019 (the green vertical line in Fig. 3). After visual inspection 
of the unbinned light curve for band 8, we found that the timescale 
of this transition is likely to be about a week or less. This timescale 
is much shorter than either the theoretical hydrodynamic times- 
cale of a few weeks to a few months” or previously observed tim- 
escales of changes in the H,O features seen in evolved red stars**’”. 
Because the gas clump cannot move outside the line of sight to the 
photosphere on such a short timescale, this transition indicates the 
occurrence of an episodic bursty event. One may speculate that the 
cause of this unusual event is related to another rare event, the Great 
Dimming; if so, it would be easy to imagine that such an unusual 
event could occur. One conceivable scenario involves a shock wave 
passing through the clump. It has been suggested that such a shock 
emerged from the bottom of the photosphere in January 2019 
and propagated from there to the outer envelope of Betelgeuse”’. 
Together with the fact that A(V) began to increase at the time of this 
transition (Fig. 2d), this shock might be related to the triggering 
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Fig. 2 | Time variation of the stellar parameters of Betelgeuse. The black 
dots with error bars represent the measured values and their s.e.m., and 
the blue lines and blue shading show versions smoothed by averaging five 
points to guide the eye and their s.e.m. a, The light curve of Betelgeuse at 
band 3, centred at 640 nm, in AB magnitudes. b-e, The results of the SED 
fitting. b, R. ¢, Toy. d, ACV). e, t19. The time of the Great Dimming is indicated 
by the red vertical line. 


process of the Great Dimming. Future theoretical investigations and 
time series of spatially resolved infrared observations would enable 
us to tackle this further mystery. 
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Fig. 3 | Molecular H,O emission and absorption. The flux excess in band 
8, where molecular H,O features exist. The blue line and blue shading 
show a smoothed version and s.e.m. by averaging five points of the original 
black dots with error bars (measured values and s.e.m.) to guide the eye. 
We found a rapid transition in the flux excess from negative (absorption) 
to positive (emission) in early April 2019 (green vertical line), which 
occurs 10 months before the Great Dimming in early February 2020 (red 
vertical line). 


Methods 


We first aimed to establish and test a new concept, using meteorological satellites as 
astronomical space telescopes, to obtain multiband, long-term and continuous light 
curves uninterrupted by the Sun. Second, we applied this technique to Betelgeuse 
to investigate its mysterious Great Dimming. 

Images of celestial objects, particularly the Moon, obtained by meteorological 
satellites have long been used only to calibrate and ensure the quality of Earth 
observations**~“’. In contrast, we attempted to use these images to investigate stellar 
astrophysics. Here, we thoroughly describe the data properties and the photometric 
methods and their validation. 


Basic properties of the Himawari-8 scanning observation. Geostationary 
meteorological satellites, such as Himawari-8, observe the entire disk of the Earth 
by sequential scans with several rows, each of which is termed a swath'*"’. Each 
swath is scanned with a one-dimensional (1D)-like detector array from west 

to east. Due to this scanning procedure, the properties of images obtained by 
Himawari-8 have several differences compared with those obtained by typical 
astronomical observations using two-dimensional devices. Nevertheless, we 

can use the reduced and distributed data, termed the L1B product, as ordinary 
astronomical images, provided that several caveats are taken into consideration. 

First, pixel navigation and data sampling are well controlled and regular, 
and their errors are usually smaller than 0.5 pix for Himawari-8". Thus, we can 
precisely transform the position in the image onto its location in the sky. 

Second, the fill factor of the pixels in the north-south direction is almost 
100%, and the east-west sampling interval is smaller than the pixel size*'. As such, 
we do not need to be concerned with undersampling due to gaps between pixels. 
Moreover, after the resampling process during the reduction described below, each 
image is transformed to the fixed square grid, which is very easy to use. 

Third, because each row of unresampled images is observed by a single 
detector, inaccurate calibration of the detector occasionally introduces a stripe 
pattern into the obtained images. The strength of such a stripe has been shown 
to be small for observations of Earth scenes”, but it has not yet been closely 
examined for faint sources such as stars. Therefore, it should be kept in mind 
that flux measurements of stars are possibly influenced by this relatively large 
systematic error. 

Finally, other possibly minor effects include quantization noise, nonlinearity 
of the detector response, and frequency of the calibration. We will further examine 
some causes of systematic and/or statistical errors and validate the photometric 
results in Supplementary Discussion. 

In this work, we limit our analyses to the images observed with the AHI 
on board the Himawari-8 satellite, which is a ‘third-generation meteorological 
satellite’ characterized by access to as many as 16 bands'’. There also are three 
other third-generation satellites equipped with 16-band imagers: GOES-16 
and GOES-17*, which are operated by the National Aeronautics and Space 


PETERS 


Administration and the National Oceanic and Atmospheric Administration, 

and GEO-KOMPSAT-2A*, which is operated by the Korea Meteorological 
Administration. Although they have been operational since November 2016, 
March 2018 and July 2019, respectively, flux values outside the Earth’s disk are 
found to be masked and unavailable in their published standard L1B products, 
which makes these data unusable for our analyses. Combined analyses with these 
satellites would improve the precision of our results, if their images of the space 
around the Earth’s disk were to be made publicly available. 

Reduction of the AHI images is performed by the Japan Meteorological 
Agency*'™’. This involves geometric regridding into a fixed coordinate system 
using a Lanczos-like resampling; a linear and global flux calibration using deep 
space as the background, with a solar diffuser for the optical and a blackbody for 
the infrared as standard flux sources; and nonlinear flux calibrations based on 
prelaunch ground tests. No corrections are made for atmospheric effects, nor is 
masking for cosmic rays performed. 


Photometry of the Himawari-8 image. About 6 yr of 16-band images have been 
obtained by the Himawari-8 satellite at a frequency of once every 10 min. We 
restricted the scope of our analyses to data observed from January 2017 to June 
2021 because systematic biases in the flux calibration and in the positioning were 
identified in observations from July 2015 to November 2016 (Supplementary 
Discussion). From this large number of images, we tried to find and measure the 
light from Betelgeuse, which is moving with the speed of the Earth’s rotation, 
15°h"!, and which is found to be located in the field of view about once every two 
days. First, we calculated the times when Betelgeuse is in the direction of the Earth 
using the Astropy package’’, and we collected all the images in which Betelgeuse 
could be detected. Then we searched for the location of the Betelgeuse signal in 
each of these images in bands 1-10, where the signal-to-noise ratio (SNR) is high 
enough to detect the star. For each of these high-SNR bands, we searched for the 
location of the star by choosing the pixel with the largest counts among those pixels 
where Betelgeuse is expected to be located. In contrast, we inferred the locations 
in the other bands (bands 11-16) by shifting from the position of Betelgeuse in 
band 7, where the SNR is the highest, by around 30 arcsec depending upon the 
differences in the times of observation for the two bands. 

The pixel scale of the published images in each band is larger than the actual 
radius of the point-spread function estimated from the radius of the telescope. 
Thus, the light from a point source occupies only 1 x 1 pix’ of the unreduced 
images, and it spreads out to ~2 x2 pix during the reduction process (Fig. la-c). 
Due to this small point-spread function radius, point-spread function photometry 
does not work well, so instead we used the aperture-photometry technique with 
a 7X7 pix’ square aperture. The background flux in the aperture was estimated 
and subtracted by fitting the flux within a square annulus around the star using 
a two-dimensional polynomial function. We set the inner width of the annulus 
to be 7 pix and the outer width to be 43, 21 and 21 pix for the bands with spatial 
resolutions of 0.5, 1.0 and 2.0km, respectively. We selected these widths to 
minimize the impact of stripe-shaped fluctuations of the background. We selected 
the order of the polynomial to be in the range 0-4 for each measurement so that 
the Akaike information criterion calculated from the fitting residuals is minimized. 
Then, we estimated the stellar flux by summing the background-subtracted 
spectral irradiances within the 7 x7 pix’ aperture. 

Since the background shapes are not guaranteed to be continuously connected 
between adjacent swaths, we excluded from our analyses cases in which the 7 x7 pix’ 
aperture spans multiple swaths. Furthermore, we trimmed the square annulus used 
for background estimation by requiring that the annulus does not span multiple 
swaths. Moreover, to avoid contamination from terrestrial light, we introduced 
distance thresholds from the edge of the Earth’s disk (Extended Data Fig. 1). 


Validation of the photometry. Here, we performed a wide variety of validation 
steps, not only with Betelgeuse but also with other bright stars located in the field 
of view of the Himawari-8 satellite (Extended Data Fig. 2). From the validation 
described in Supplementary Discussion and Supplementary Figs. 1-6, we found no 
evidence of any serious systematic bias in our photometry, and we thus conclude 
that the measurements are ready for further analyses. Nevertheless, it would be 
important, in a coming decade, to examine whether or not there are overlooked 
systematic effects in our method. 


SED fitting to the Himawari-8 light curves. Figure 1d,e shows the SED averaged 

in time together with a model SED having the typical atmospheric parameters 

of Betelgeuse. The photometric bands with wavelengths of <8 j1m (bands 1-10) 

are almost unaffected by the circumstellar dust emission, and thus they can be 

used to determine the atmospheric parameters of Betelgeuse. In contrast, the 

other bands receive radiation from silicate-dust emission at around 101m, and 

they can thus be used to measure the amount of dust formed around Betelgeuse. 

We therefore first fitted the SED for bands 1-10 with a model photosphere for 

Betelgeuse to determine the photospheric parameters—R, T, and A(V)—and we 

then determined 7,) using bands 12-14, which have the highest SNR for the dust 

emission. Five examples of the fitting results are shown in Extended Data Fig. 3b-f. 
Our final time-series photometric catalogue of Betelgeuse is an incomplete one; 

the catalogue lacks fluxes in several bands on several observations. 
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This is especially serious for longer-wavelength bands, where the sky-background 
emission is stronger, and thus a larger number of photometric measurements was 
removed from the catalogue when the star was located near the edge of the Earth’s 
disk. Moreover, the signal-to-noise ratios in the mid-infrared bands are typically 
small, ~1. To deal with these difficulties, we binned each light curve over 5-12 
observations, requiring that all the binned data points for bands 1-10 include 

five or more points. Note that the number of days in each bin is automatically 
determined by the observation frequency, which depends on the observation 
epoch. For example, there is a clump of observations in May-June 2018 with a 
frequency of 0.85 d, which leads to the smaller number of days in the bins. Finally, 
we calculated the resulting median and its standard error for each bin. 


Model SED for the photosphere. We fitted each point of this ‘complete’ SED curve for 
bands 1-10 with a photospheric model SED using the LMFIT package*. To obtain 
the model spectrum for the photosphere, we interpolated a grid of BT-NextGen 
model spectra” with solar metallicity and the surface gravity log g=0.0, assuming 
that this 1D model can be used to model the spectra of Betelgeuse in any phase. We 
used Cardelli’s reddening law® with the selective-to-total extinction ratio Ry=3.1 
(ref. *'), assuming that this law can be extrapolated to longer wavelengths. We also 
used the parallax of 6.55 mas measured with the Hipparcos satellite’. Although 

the accuracy of the parallax of Betelgeuse is controversial**’, the choice of parallax 
value only changes the fitted radius in a relative way and does not affect the 
interpretation of the results. For the same reason, we do not take into account the 
error in the parallax value. The model spectra are then convolved with the response 
function of the AHI (September 2013 version retrieved from https://www.data. 
jma.go.jp/mscweb/en/himawari89/space_segment/spsg_ahi.html) to calculate the 
model SEDs for Betelgeuse, and we fitted the observed SED with the model. 


Dependence on the model selection. For the grid of model spectra, we first tried 

the MARCS grid® with 15M, (refs. **°**°), which is usually used to model the 
spectra of red supergiants’'’. We found that the MARCS grid gives a smaller 
Akaike information criterion than the BT-NextGen grid, probably because the 
wavelengths and the strengths of the molecular bands in the MARCS grid are more 
accurate. Moreover, considering the fact that several previous works”! succeeded in 
reproducing observed optical spectra of red supergiants dominated by molecular 
TiO bands with MARCS models, MARCS models appear to be an appropriate 

grid for the synthesis of TiO bands in red supergiants. However, we found that the 
results of the fits obtained using the MARCS grid were physically unreasonable 

for most cases that we considered, for example, A(V) <0 mag, although the reason, 
possibly related to 3D effects and/or circumstellar gas, is unclear. Thus, instead 

of the MARCS grid, we used the BT-NextGen grid to compute the model SED. 
Testing both the MARCS and NextGen grids, we found that the model parameters 
we obtained using the BT-NextGen grid vary with time in the same way as 

those obtained using the MARCS grid; that is, the differences between the time 
variations obtained using the two grids are roughly constant (Extended Data Fig. 
4a-c). Therefore, the results of this paper using the NextGen grid are not affected 
quantitatively by the choice of model grid, as long as the time variations are 
interpreted in a relative way. 

We have also tested whether or not other model assumptions affect the results. 
For example, a change in log g from 0.0 to —0.5 makes T., systematically ~80 K 
warmer and A(V) systematically ~0.3 mag larger. Similarly, changing Ry from 3.1 
to 4.4 (ref. °') decreases T, by $10 K and A(V) by $0.05 mag. From these tests, 
we concluded that, although the absolute values of R, T,;, and A(V) depend on the 
model assumptions, their relative variations appear to have been determined to an 
accuracy of ~10K in T,, and ~0.05 dex in A(V), both of which are much smaller 
than the fitting errors, ~40 K and ~0.2 dex, respectively. 


3D effects in the photosphere. The temperatures that previous works*'*°! and we 
have both derived assumed a 1D geometry for the model atmosphere. Thus these 
temperatures are affected by 3D effects in the optical spectrum of the real star’. 

In fact, the SED of a model with a given T,;, especially the strength of the TiO 
molecular absorption bands that appear in the optical spectra of red supergiants', 
has a variety of shapes due to 3D effects*. Therefore, when we fit the observed SED 
of the real 3D photosphere of Betelgeuse with a 1D model SED, the fitted values 

of T,; and A(V) depend not only on T,; and A(V) themselves but also on the 3D 
effects*’. In particular, a time variation of the patchiness of the surface” could 
affect the time variations of T,,,and A(V) (ref. ’). 

As a simplified test of the 3D effect on our fitting result, inspired by the optical 
image of Betelgeuse obtained during the Dimming'”’, we considered a model 
photosphere composed of two regions with different T,;,, which is similar to the 
model used in a previous work’. In this model, we assume that the cooler part of 
the photosphere has a varying temperature and a varying area and that the warmer 
part has a fixed temperature of 3,650 K. In contrast to the previous model’, we are 
able to consider two additional free parameters, A(V) and R; without these two 
parameters we could not obtain a satisfactory fit to our multiwavelength data. 
Altogether, we considered four free parameters (R, A(V) and the temperature 
and area fraction of the cooler part); for example, when the entire area of 
the photosphere is occupied by the cooler part, the fraction is 1.0. With this 
four-parameter model, due to the strong degeneracy between A(V) and the area, 
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sometimes the fitting failed and gave unphysical results. To remove these cases, we 
excluded fitted results with areas greater than 1.0 or less than 0.0 as well as results 
for which the temperatures reached the boundaries of the acceptable range (we 
assumed 2,800 and 4,000 K for the lower and upper boundaries, respectively). 

Extended Data Fig. 4d-g compare the fitted results obtained with the 
three-parameter model adopted in the main text and the four-parameter model 
described here. We found that with the four-parameter model the contribution of 
the increased A(V) to the Dimming decreases but that it is still larger than zero 
(~0.4mag). This result supports the robustness of our results in a qualitative sense. 
Since the four-parameter model examined here is too simple to take 3D effects 
fully into account’, future theoretical modelling with self-consistent hydrodynamic 
models” will be necessary. 


Comparison with previous radial-velocity (RV) measurements. To check our SED 
fitting results, we compared the radius determined herein with previous RV 
measurements. Because the radius is given by the integration of the RV, we can 
validate the radius curve using literature RV curves. 

First, we found time lags of up to ~80d between the times of maximum 
brightness and minimum radius (Fig. 2). Although such a large lag was not 
expected in a simple 1D hydrodynamic model”, a lag between light and RV 
curves has been observationally found for Betelgeuse”’»°°. Moreover, we found 
clockwise trajectories between the R/T, and band 3 flux for most observation 
periods (Extended Data Fig. 3h—-o), which has been predicted theoretically™”’. The 
exception, the relation between R and band 3 flux between 27 October 2017 and 
24 December 2018 shown in Extended Data Fig. 3i, could be attributed to effects 
other than the fundamental-mode pulsation, such as the long secondary period”; 
however, the actual cause remains unclear. 

Second, we estimated the RV curve by calculating the differential of our radius 
curve, and compared it with literature RV curves'*"°. This comparison reveals an 
overall good agreement between our data and the results cited above, supporting 
the validity of our radius measurement, and hence other parameters. However, 
minor differences of up to ~5kms™ may exist, particularly in September 2018- 
November 2018 and September 2019-December 2019. Plausible causes of these 
differences are (1) our imperfect modelling of Himawari-8’s SEDs with a 1D model 
atmosphere, as discussed in the previous subsection, and (2) differences in the 
atmospheric layers that were observed in the present and previous studies, and 
hence the real, physical difference in RVs”. 


Model SED for the dust emission. After fitting the SED curve for bands 1-10, we 
calculated the expected photospheric flux curves for bands 11-16 and subtracted 
them from the observed flux curves. These residual fluxes are likely to originate 
from the circumstellar dust emission, and they are useful for estimating the 
variation of the amount of circumstellar dust. 

For the model spectra of the circumstellar dust emission, we used the 
radiative-transfer code DUSTY V2”. For the settings of the DUSTY calculation, we 
basically followed the procedure described in a previous work”, except for several 
modifications. First, we only fitted z,, and fixed the other parameters. We used the 
‘warm silicates’ dust composition”', a constant grain size of 0.3 pm and a constant 
dust temperature of 1,300 K at the inner boundary, which roughly corresponds 
to the sublimation temperature of silicate dust**. Although inner-boundary 
temperatures less than 1,300 K are sometimes found using observations of red 
supergiants”, we adopted this value since a different inner-boundary temperature 
only results in a systematic bias in 7,, as long as the temperature is independent 
of time. We also tried the assumption of a constant radius at the inner boundary, 
and we found that the results do not change in a qualitative sense. For the stellar 
radiation that illuminates the dust, we used the time variation of T,, determined 
in the previous section to calculate the natal photospheric SED, after applying 
spline smoothing to reduce the noise in each measurement. Then, we calculated 
the expected silicate-dust emission in bands 12-14 on the 7,, grid and determined 
the 7,, value that best matches the observed flux-residual curve for each of 
bands 12-14. Finally, we calculated the mean of the three 7, values determined 
from bands 12-14 to obtain the final estimate of 7,,. We note that the z,, that we 
determined is in fact not the actual optical depth at 10 jm; it is rather a model 
parameter that we use to match the observed and model fluxes of bands 12-14, and 
we use it as an indicator of the amount of circumstellar silicate around Betelgeuse. 


Data availability 

Reduced images taken by the Himawari-8 satellite are publicly available at the Data 
Integration and Analysis System via https://diasjp.net/en/. The final photometric 
product of this paper is available at https://d-taniguchi-astro.github.io/homepage/ 
Data_Himawari_en.html, which is archived as Supplementary Data. 


Code availability 

We used the LMFIT package* to fit the observed SED in the optical to near 
infrared range with publicly available NextGen* and MARCS* model grids. We 
also used the DUSTY V2 code” to calculate the dust emission. We have not made 
the codes for the photometry publicly available because they are not prepared for 
open use. Instead, we have made the photometric catalogue available, as described 
in the data availability statement. 
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Band Xo" AX Pixel scale Man?  Finaxe LSB° Zmint LB test p value® 

No. (yum) (um) (km)? (arcsec) (mag) (nW m~? pm) (km) Rigel Procyon 
1 0.471 0.041 1.0 5.76 2.4 597 0.38 75 0.419 0.509 
2, 0.510 0.031 = 1.0 5.76 2.5 560 0.35 75 0.174 0.210 
3 0.639 0.082 0.5 2.88 2.4 121 0.31 75 0.003 0.321 
4 0.857 0.034 =1.0 5.76 0.8 288 0.182 60 0.884 0.859 
5 1.610 0.041 2.0 11.53 0.7 287 0.045 65 0.341 0.204 
6 2.257 0.044 92.0 11.53 0.9 89.1 0.0141 60 0.938 0.609 
7 3.89 0.20 2.0 11.53 2.0 40.2 0.000 79 55 0.964 0.329 
8 6.24 0.82 2.0 11.53 0.3 29.9 0.0048 90 0.331 0.836 
9 6.94 0.40 2.0 11.53 —0.1 36.4 0.0058 100 0.288 0.396 
10 7.35 0.19 2.0 11.53 —0.3 38.8 0.0031 100 0.893 0.382 
11 8.59 0.37 2.0 11.53 —0.7 50.6 0.0040 100 0.270 0.500 
12 9.64 0.38 2.0 11.53 —0.8 49.7 0.0039 130 0.689 0.316 
13 10.41 0.42 2.0 11.53 —1.2 47.5 0.0038 90 0.756 0.870 
14 11.24 0.67 2.0 11.53 —1.0 44.3 0.0035 80 0.913 0.250 
15 12.38 0.97 2.0 11.53 —1.4 39.4 0.0031 100 0.658 0.285 
16 13.28 0.56 2.0 11.53 —2.5 36.0 0.0058 100 0.351 0.074 


Extended Data Fig. 1| Basic properties of the Himawari-8 AHI images and the thresholds used for photometry. @Each band is centered on A, and has the 
FWHM Aj; Spatial resolution at the sub-satellite point, i.e., 35,786 km from the satellite; “Sky-3o limiting AB magnitude for each band at the declination 
of Betelgeuse; (Maximum valid flux per pixel in the calibrated L1B file format of Himawari-8 images. According to results of ground tests”, the saturation 
flux of the AHI itself is larger than these values; °Least significant bit (LSB) of L1B images; i.e. the flux of each pixel is measured in increments of LSB; 
‘Minimum distance between the Betelgeuse line of sight and the edge of the Earth's edge accepted for photometric analyses in this study; 8p value of 
Ljung-Box test’, with the time delays up to 30 days for the two non-variable stars, Rigel and Procyon. 
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Name HD Sp. Type* Dec.* Frequency” 
number (deg)  (days~+) 
Betelgeuse 39801 MI1-M2Ia-Iab +7.41 0.58 
Rigel 34085  B8lae —8.20 0.91 
Procyon 61421 FSIV-V+DQZ +5.22 0.33 
6 Oph 146051 MO.SHI —3.69 0.29 
a Cet 18884 M1.5Iia +4.09 0.48 





Extended Data Fig. 2 | The stars for which we performed photometry. *Spectral type and declination (J2000.0) taken from the SIMBAD database” on 15 
June 2021; *Typical frequency at which the star appear in the Himawari-8 images. 
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Extended Data Fig. 3 | Results of the photometry and the SED fitting for Betelgeuse. a, Light curves of Betelgeuse in all bands. The black dots represent 
individual measurements, and the blue circles and vertical bars represent the medians and standard errors after binning, respectively. The photometric 
measurements in the gray-shaded areas, that is, in 2015 and 2016, were not used in the analysis. b-f, Five results of the SED fitting as examples, whose 
observation epoch are shown in g. The red and brown circles represent the observed and fitted spectral irradiances for each band, respectively. The 
purple and blue lines represent the fitted spectrum before the convolution with the response function. g, The light curve of Betelgeuse at Band 3 in AB 
magnitudes, as in Fig. 2. h-k and I-o, The relation between the radius R or effective temperature T,,,, respectively, and Band 3 magnitude with the smoothed 


trajectory created by averaging five points. 
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Extended Data Fig. 4 | Comparison of the results of SED fitting with different model assumptions. In all the panels, the black circles indicate the 
results adopted in this paper (the same as Fig. 2), which use the grid of Next-Gen model spectra*’, assuming three free parameters in the fitting. a-c, The 
effect of the adopted model-spectrum grid: Next-Gen”? (blue circles) vs. MARCS*® (upward-pointing orange triangles). d-g, The effect of the number of 
parameters: three parameters (blue circles) vs. four parameters (downward-pointing pink triangles). 
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